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FOREWORD 


This, Volume I of the report "High Speed Cylindrical Roller 
Bearing Analysis," documents the analysis, program design and 
algorithm detail employed in the generation of the computerized 
analytic design tool CYBEAN . All efforts involved in the gene- 
ration of the code were sponsored by the NASA-Lewis Research 
Center of Cleveland, Ohio, under the administration of the 
Fluid System Components Division. The technical monitor was 
Mr. H. Coe. The work was performed under Contract No. NAS3-20068 
at SKF Industries, Inc., King of Prussia, Pennsylvania, during 
the period September, 1976 through July, 1978. 

Technical project leadership was executed by Mr. R. J. 
Kleckner, with contributions from: Drs . V. Castelli and 

J. Pirvics, Messrs. W. J. Crecelius and M. Ragen and 
Ms. M. M. Din on. 


VOLUME I 


AL78P022 


HIGH SPEED CYLINDRICAL ROLLER BEARING ANALYSIS 
I. INTRODUCTION 

Engine performance which satisfies the anticipated effi- 
ciency, speed and specific thrust requirements for aircraft in 
the coming decade project increased demands on their component 
performance. Engines will be required to operate at higher 
temperatures and speeds and to do so while maintaining low over- 
all mass and high structural integrity. Bearings in particular 
emerge as critical system design elements. 

Temperatures, which will be increased to extract more of the 
available thermodynamic engine cycle efficiency, affect bearing 
performance by threatening the dimensional stability of the 
system. Thermally induced excursions in geometry must be com- 
pensated for within the original bearing design to avoid pre- 
mature and possibly catastrophic system failure. 

Engine shaft speeds will be increased to derive greater 
compressor efficiency. Increased speeds, however, accentuate 
centrifugal effects. Dominance in raceway loading is trans- 
formed from inner to outer rings. Failure to maintain loaded 
contact on the inner ring across th& spectrum of operating 
speeds results in the increased hazard of roller skidding, 
unnecessary heat generation and unstable performance. Combined 
with increased shaft diameters, required by system stiffness, 
increased speed raises questions beyond the limits of currently 
available bearing design analysis. This is especially true for 
cylindrical rolling element bearings where the diameter (mm) _ 
speed (rpm) product, DN, approaches 3.0 x 10^. 
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Field experience supports the assessment that available 
understanding and design tools are inadequate. Roller deterio- 
ration in bearings operating in the noted DN regime provides 
the physical evidence. This deterioration indicates that roller 
skidding occurs in lightly loaded high speed applications. 

High bearing temperatures and erratic kinematic behavior 
contribute to bearing degradation. Further examination 
indicates that reliable bearing operation depends on the ability 
to quantitatively assess the extent to which rollers will skew 
and which of the available design parameters can be manipulated 
to control this particular roller orientation during operation. 


The deficiency in current design analysis tools is that 
the anticipated engine performance demands require the under- 
standing of behavior which has previously resided in the realm 
of secondary effects. 


This emergence of effects from secondary to primary sig- 
nificance has dictated the construction of new analysis software 
The impetus for such effort has' been further accelerated by the 
high labor and material costs required for purely experimental 
hardware research. 
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This report details the work performed to permit computation 
of the behavior of cylindrical rolling element bearings under high 
speed and misaligned shaft conditions. The program objectives were 
pos'ed to enable accurate roller skew prediction under a wide 
variety of design geometry and loading conditions. Emphasis 
was also placed on the accurate assessment of general geometric 
preload design. Finally, it was recognized that such detailed 
information was deceptive without an integral and system related 
thermal analysis. 

The work statement contained in Appendix A details the 
specific objectives relating to this NASA sponsored project. 

The material which follows details the procedures used in 
attaining the stated goals and instructs the reader in the use 
of the generated computer program CYBEAN. The documentation 
is presented in two volumes. Volume I describes program archi- 
tecture and follows with modular documentation detailing 
formulation and mathematics. Volume II is written as a User's 
Manual for CYBEAN. The reader is presented with material which 
guides sensible and accurate use of the program. Volume II 
appropriately references Volume I for algorithm background or 
the details of simulation mathematics. 

II. PROGRAM ARCHITECTURE 

CYBEAN has been structured to be a bearing analysis tool 
by the coordinated execution of modules which perform specific 
analytic tasks. 

The basic structure of the program is seen in Figure 1. 
Computation economies are accomplished by placing constants in 
structure locations external to repetitive operations. 
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CYBEAN initiates computation with geometry, lubricant and 
material definitions. The first repetitive loop is entered 
immediately thereafter, addressing the computation of interactive 
geometry fits. Arbitrary housing and outer ring profile inter- 
actions are considered here. Load distribution and displacement 
values are obtained for the specific mode of bearing support 
specified by the user. The evaluation of the field equation 
set for a bearing which contains n rollers is now started. 


The set of values for the independent variables which at 
a particular instant in time simultaneously satisfy all 
equations of equilibrium are considered to be a solution. 


These values are obtained by solving as many as 6n r + 5 simul- 
taneous nonlinear equations using a Newton- Raphson algorithm. 


The individual positions of the rolling elements within 
the complement at that instant are defined by^reference to 
the orientation azimuth angle, ipi, of the first rolling 
element. 


Because of the nature of the governing equation set, an 
accurate initial guess of the variable values would result in 
considerable savings in computer execution time. However, the 
nonlinear coupling of a large set of unknowns within the phy- 
sical system creates difficulty in the generation of an accurate 
initial guess. 

To overcome this problem and maintain computation time 
at an economic level, the simpler elastic problem is solved 
first. The set of variable values which satisfy the equilibrium 
of elastic forces are then used as an initial guess; for the 
more complicated solution; which includes friction. This 
solution is also used to determine the equilibrium load 
distribution which defines housing fit pressures. The basic 
problem is described by Liu [1], the number of degrees of 
freedom consider d are 2n r +2. 

. . 4 
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The general field equation set is partitioned prior to 
its evaluation and eventual solution. Independent variables 
describing the inner ring position and roller radial position 
are addressed in conjunction with field equations governing 
ring and roller radial equilibrium. Equation cross coupling 
is established through the roller radial position. This subset 
is brought to equilibrium by Newton- Raphson iteration. 

Once inner ring equilibrium. is satisfied, the field equa- 
tion set is again partitioned to include roller axial, radial, 
skew and tilt moment equations of equilibrium. Coupling is 
maintained through the associated displacement variables. At 
this portion of the analysis, roller end flange contact forces 
are equilibrated with raceway induced loading. 

A final partition is made such that the. field equations 
defining force equilibrium in the z-direction, moments about 
the roller axis, radial and cage equilibrium are considered. 
Associated variables are: radial position, rotational and orbital 
speeds and cage displacements. That set of independent variables 
whose values simultaneously satisfy these equations is considered 
a solution.: 

Satisfactory convergence of the iterative procedure allows 
the definition of results for a specific temperature. The 
remaining program structure addresses the computation of steady 
state and the time transient thermal performance of the 
system. 
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FIGURE 1: BASIC PROGRAM STRUCTURE 
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FIGURE 1 (CONT'D): BASIC PROGRAM STRUCTURE 
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III. ANALYSIS 

The preceding de$cription of program architecture serves to 
provide an overview of CYBEAN structure. The analysis repre- 
sented by this computation algorithm is described in the material 
which follows. 

r . . : t ' 

The bulk of the calculations performed address repetitive 
computation of displacements and forces resulting from solid- 
solid and solid- fluid interactions. These interactions are 
permitted to occur in a thermally varying environment and are 
viewed as part of a system in quasi-static thermodynamic equi- 
librium. The system is defined by geometric and material de- 
scription of the bearing complement and the rheology of the 
..lubricant used. It is considered to be driven by an external 
forcing function which consists of a constant vector set repre- 
senting definitions of shaft rotation, speed and applied radial 
force. Additional system stimulation is accepted from heat 
generated within the boundaries of the defined system volume. 

The analysis consists of an ordered series of mathematical 
descriptions which detail the piecewise and ultimately continuous 
equilibrium achieved by the system in response to the applied 
forcing function. It describes the response of the system in 
terms of force, displacement and temperature distributions. 

The bearing analysis is initiated through the definition 
of the coordinate space in which it is described. An inertial 
reference frame, chosen to be coincident with an inner ring 
reference frame when the bearing is unloaded, is established 
such that its X-axis defines the bearing longitudinal centerline. 
Each element is then assigned a moving coordinate frame. For the 
rollers, this frame is coincident with an analogous inertial 
frame when the bearing is unloaded. Thus defined, a general 
displacement (composed of both translation and rotation) of 
any bearing component can be written in the terms of a coordinate 
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transformation from its inertial reference frame. Six coordinates 
for a given transformation are required to establish the displaced 
positions of each bearing component. However, if axial symmetry 
exists constraining assumptions can be imposed on the nature 
of the displacement, and the number of degrees of freedom required 
to define the rigid body displacements can be reduced. 
Specifically: 

Inner Ring - The study of misalignment effects on bearing 
performance is made convenient by the ability to define 
the misalignment angle. An additional constraint on the 
axial displacement and the presence of rotational symmetry 
can reduce the degrees of freedom to two. Thus, a very 
restricted representation describes translation of the 
ring within its radial plane, 

Cage - Motion can be restricted as with the inner ring, 
although rotation is permitted about its longitudinal axis. 

Rollers - Axial symmetry and the recognition that trans- 
lation in the direction tangent to .the pitch circle is 
negligibly small can limit the motion to 4 degrees of 
freedom. 

Those quantities which define the degrees of freedom through 
which the coordinate frames may move form a subset of the inde- 
pendent variables in the field equation set. The complete set 
includes the roller rotational and orbital speeds. 

The definition of general coordinate space transformations 
and the general geometric description of each element in the 
bearing complement (within its particular coordinate frame) 
eases the formulation and permits flexibility in the numerical 
experimentation required to attain optimum solution algorithms. 
Appendix B details this basic transform structure. 
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The effort has been directed at formulating the description 
to allow full generality. This modular freedom is guaranteed 
by the automated coupling provided by the transformations noted 
above . 

The geometric characterization of rings, rollers and cage 
is elaborated upon in Volume II of this work. Some of the figures 
used in Vo . ume II are repeated here, in Appendix C. 

To introduce physical significance, spacial definition 
of the system and its components requires characterization 
by specific material properties and lubricant rheology. 

Solids are characterized as linearly elastic deforming 
materials whose properties are temperature independent. Spe- 
cific material properties 1 are found in Volume II. 

Lubricating fluid is characterized as having properties 
that vary with both temperature and pressure. Models used to 
evaluate properties are shown in Appendix D. 

Once geometry and material properties are defined specific 
interaction mechanisms must be described to formulate the force 
equilibrium field equation set. Two major descriptive areas 
arise: solid-solid and solid- fluid interactions. 

The solid-solid interaction is characterized as follows: 
given a set of values for the independent variables, positions of 
all the bearing components are established in space. With these 
positions fixed, the relative positions of the bearing components 
are examined to determine if the bodies contact or a free space 
exists between them. Locations and magnitudes of interpenatration 
are recorded. These values are subsequently used to determine 
contact stress, contact size and the associated elastic force. 
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The solid- fluid interaction occurs in four ways: 
raceway /rolling element traction, flange/rolling element traction, 
roller drag and cage related friction. Again, with the positions j 

i 

and speeds fixed, these effects can be computed for the raceway j 
usjing established models. In this work, models used to treat 1 
flange/ rolling element and raceway /rolling element interactions 
are identical. Roller drag is approximated as a function of 
lubricant properties, speed and geometry, using established 
equations for external flow. These interactions and the cage 
model are detailed in Appendix E. 

The analysis of induced preload through variable geometry 
ring interactions with a variable geometry housing is described 
in Appendix F. 

The completion of the description detailing interaction mecha- 
nisms allows the generation of the quasidynamic force equilibrium 
equations and attendant strategies for solution. A Newton- Raphson 
iterative procedure noted in the preceding section is used. 
Modification of the equation set through the incorporation of 
kinematic equations is detailed in Appendix G. 

The culmination of the iteration procedure which satisfies 
the j quas idynamic equilibrium field equation set permits the 
display of values for the vector of unknowns and the computation 
of bearing fatigue life. The latter is computed according to 
standard Lundberg-Palmgren methods detailed in Appendix H. 

The description to this point has addressed the bearing 
analysis performed at a given set of component temperatures. 
Thermal effects can be computed in either of two modes. 
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The first mode details the steady state operating temperatures 
of the bearing and its environment. CYBEAN will formulate the 
conservation of energy equations for ail equivalent nodal model 
of the bearing and surrounding hardware. This model is prepared 
by the user according to methods shown in Volume II. Equations 
are solved for the steady state temperature distribution. 

The second mode details time transient effects. Here, CYBEAN 
solves a system of first order nonlinear differential equations. 
Typically, initial values are taken as the solution to the steady 
state analysis. Modelling assumptions are presented in Appendix 
I, and also Volume II of this work. 

IV. CONCLUSIONS AND RECOMMENDATIONS 

The material; presented in the preceding chapters documents 
the creation of the state of the art high speed cylindrical 
rolling element bearing analysis /design tool CYBEAN. The for-, 
mulations are general and the program architecture modular to 
permit easy maintenance and expansion. 

The current edition of CYBEAN enables detailed examination 
of bearing performance and in particular permits exploration 
of causes for and consequences of bearing skew. The very gen- 
eral capability for calculation of geometrically induced pre- 
load allows design assessment of bearings for engine mainshaft 
support. The integral thermal analysis capability obsoletes 
previous isothermal computation algorithms for state of the 
art 3 x 10^ DN bearing application investigations. 
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During the course of the work performed, a perspective 
was gained on the additional requirements for the accurate simu- 
lation of design components for future generation engines. These 
requirements relate to the previously noted emergence of physical 

•' i 

phenomena from secondary to primary effect significance as oper- 
ating conditions specify higher speeds and temperatures. Proper 
assessment of influence on bearing performance demands an invest- 
ment in basic analytic formulation and corresponding detailed 
experimental verification. It is specifically recommended that: 

1. The EHD model be reformulated. Experiments leading 
to the current version, developed in the early 1970's, 
excluded regimes which are forecast for the 1980 's. 

- This inadequacy is experienced when performance 
includes a high degree of skidding. 

2. The effects of lubricant distribution be investigated. 
Thermal/dimensional stability and thus power loss in 

a bearing is exceptionally sensitive to the assumptions 
concerning the lubricant distribution within the 
bearing cavity. This sensitivity is accentuated with 
increased speed and temperature. Both of these para- 
meters have increasing effects with bearing operation 
at elevated DN values. 

3. Surface traction effects be related directly to 
material failure and thus bearing failure predictions. 

4. The microscale topography of surfaces be quantified 
in three dimensions to enable rational formulation 
•of upgraded EHD models which address physical events 
in this microscale realm. Current procedures are 
limited by excessive measurement times and mechanical 
filtering within the frequencies of interest. 
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5. The mathematics of numerical methods be more 
fully expolited to enable economically justifiable 
computation in parametric explorations. 

6. Extensive work be performed to obtain a quanti- 
tative understanding of phenomena associated with 
partial film and lost lubrication operating conditions 

The current edition of CYBEAN as well as the state of the 
art itself would benefit from investigations in any and all 
of the six items noted above. In addition, specific immediate 
modifications and capability expansions for CYBEAN warrant 
attention. These are represented as the following tasks: 

1. Expand CYBEAN to accept external shaft applied 
axial load. 

2. Expand CYBEAN to couple with a shaft and another 
rolling element bearing. 

3. Expand the fit routines to accomodate hoop 
strain. 

.• i V- t t ' ' : ' ' - «-c*‘ '/ ' f I f !».*. f'-- ' • ■ ' , , . -'.'V-- V- 

4. Modify CYBEAN to enable temperature dependent 
material property definition. 

5. Modify CYBEAN to automate the partitioning 

algorithm to respond to operating vector 
extremes . ' i 

t 

6. Expand CYBEAN to contain line plotter printing 
capability for time transient run information 
display. 

7. Expand CYBEAN to include a generalized interative 
thermoelastic outer ring housing description. 
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, In conclusion a recommendation is made for the expanded 
parametric exercise and documentation of obtained results. 

These results and the current edition of CYBEAN should be incor 
porated within a design guide for cylindrical rolling elements 
to provide a practical tool for the designer and analyst of 
reliable aircraft engine bearings operating in and beyond the 
3 x IQ** DN regime. 
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The Contractor shall furnish the necessary personnel, materials, 
services, facilities, and otherwise do all things necessary 
for or incident to the completion of the work set forth below: 

Task I - Program Development 

1 

I 

The Contractor shall design, develop and furnish a fully opera- 
tional computer program that will predict the performance charac- 
teristics of high-speed cylindrical roller bearings, with the 
following requirements: 

A. Allithe capabilities present in the 

bearing subprogram contained in SKF computer program 
No. AT76Y001 currently installed at NASA-Lewis Research 
Center shall be contained in this new program. 

B. The new program shall contain a "Roller Skew Analysis" 
wherein each roller shall be permitted to skew and 
this skew shall be resisted by flanges on either the 
inner or outer ring. 

C. Allowable program input shall include distinct roller- 
end spherical radius, flange angle, flange profile 
radius and roller- flange axial play. 

D. Each roller shall per permitted six degrees of freedom- 
these are : 

1. Radial displacement 

2. Axial displacement 

3. Angular displacement about the roller Z axis 
(Roller Tilt) 

4. Angular displacement about the roller Y axis 
(Roller Skew) 

5. Rotational speed (about roller X-axis) 

6. Orbital speed (about bearing center) 

E. The bearing shall be permitted to have a complement of 
up to 50 rollers . 

F. The computation of thq set of values of the variables 
in I - D that results in roller quasidynamic equilibrium 
shall be calculated through consideration of the 
elastic and friction forces imparted to each roller by 
the raceways, free lubricant, cage, and by the ring 
flanges. Inertial loads shall also be considered. 
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G. Important parameters such as EHD film thickness and 
Hertz stress shall be determined for each roller 
end- flange contact. 

H. Inasmuch as certain computations will be similar to 
previous programs, existing subroutines shall he used 
as much as possible. In particular, the following 
SKF subroutines shall be available for use with the 
new program: 

PROPST , LUPROP, LUBCON, GUESCG, FIT, INTFIT , SI'M'EQ, 
SOLVXX, INSOLV, PARDER, EQCHEK, EQWRIT, ERCHEK, LIFE, 
LRHS , BOT, VISCOZ, ALPHAO, DRAGNO, STCON, STARFC, 
HOHI, HDFRIC , HRFRIC, EHDSKF , FRICTN, NASEHD , DRAG, 
CAGESP, CAGEQ , C GLAND, CGWET, CGRE, CGBALI, CGNRMN, 
CGEHDD, CGHDP , CGFRN , CGHDF, CGEHF, SUMCGL, SUMRE, 
BRCAX, BRAX, TEMPIN, INPUM, RWHTC , RWGRWHC, TMAP, 

NET, NETEET and STEPLMA. 


I 


The computer program shall be capable of treating each 
of the following four methods of using flexible out- 
of- round bearing ring designs to induce roller preload 



The outer 
able wall 


ring has an elliptical OD and a vari- 
thickness (round ID). 


2. The outer ring- has an elliptical OD and a uniform 
wall thickness (elliptical ID) . 


3. The outer ring has a round OD and a variable wall 
thickness : (elliptical ID) . 

4. The outer ring has a three-lobed OD and ID with 
a uniform wall thickness. 


J. The new program shall be designed for eventual instal- 
lation into the SHABERTH/NASA SKF Shaft-Bearing 
Computer Program No. AT76Y001 as a subprogram module. 

K. The program shall include a capability to predict 
steady state and transient thermal performance of the 
roller bearing system. 

1. The mechanical performance analyses shall inter- 
act with the temperature mapping-heat dissipation 
analyses so that 

(a) The computed bearing component temperatures 
will influence bearing diametral clearances 
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(b) The computed lubricant temperatures will 
influence lubiicant viscosity, density and 
pressure -viscosity coefficient and 

(c) The calculated bearing frictional heat 
generation will influence the thermal 
analysis. 

; I . . ; 

2. The thermal analysis shall include the capability 
for calculation of roller- raceway and roller- 
flange heat transfer coefficients. 

Task II - Program Demonstration 

The Contractor shall demonstrate the use of the program developed 
in Task I by performing a computer analysis of a cylindrical 
roller bearing having specifications and operating conditions 
as described in Attachment 1. 

A. The purpose of the first study shall be primarily to 
demonstrate the ability of the program to predict 
roller skew and the resulting roller-end, flange inter- 
actions. Therefore, to accomplish this, two sets of 
roller end, flange geometry shall be examined under a 
condition of aligned and misaligned bearing rings, as 
noted in Attachment 1. The thermal analysis shall not 
be used in this study. 

B. The purpose of the second study shall be primarily to 
demonstrate the ability of the program to predict 
roller bearing performance as a function of the cal- 
culated temperatures. The steady state thermal analy- 
sis shall be used. To accomplish this, one set of 
roller end, flange geometry shall be examined under 
one condition of aligned and one condition of mis- 
aligned bearing rings as noted in Attachment 1. 

Task III - Program Installation 

The computer program shall be installed on the UNIVAC 1100/40 
computer at the NASA Lewis Research Center. 

A. The Contractor shall supply a program user's manual 
which will demonstrate and review program input, output, 
structure and solution methods. 

B . The Contractor shall supply one copy of the program 
source listing. 
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C. The program's satisfactory installation at Lewis t^ill 
be confirmed when the source program supplied by the 
Contractor has been successfully compiled, and executed, 
and the output correlated with the results of a test 
case (one of the data points from Task II) provided 

by the Contractor. 

D. As part of the requirement of this contract, the 
Contractor shall provide qualified technical assistance 
for a period not to exceed forty (40) engineering man- 
hours, to advise the Lewis Project Manager of any 
deficiency which may preclude proper operation, in the 
NASA- Lewis Research Center computer, of the computer 
program to be delivered hereunder. 

Task IV - Reporting Requirements 

Technical, schedular, and financial reporting shall be in accor- 
dance with the Reports of Work attachment (with NASA Forms 533M 
and 533P) which is hereby made a part of this contract except 
as modified below: 

A. The number of report copies to be submitted is as 
follows : 

1. Monthly Technical Report, a maximum of twelve (12) 
copies. 

2. Contractor's Financial Management Performance 
Analysis Report (533P) (Blocks 1, 2, 6 and 11 only), 
maximum of twelve (12) copies. 

3. Monthly Financial Management Report (NASA Form 
533M) , a maximum of four (4) copies. Column 8(a) 
slhra.ll contain the cost estimates for the month 
following that reported in column 7(a). Column 
8(b) shall contain the cost estimates for the 
month following that reported in column 8(a). 

B. The reporting categories to be reported on the 
Contractor's Monthly reports are as follows: 

Reporting Category Description 


Task 

I 

Program 

Development 

Task 

II 

Program 

Demonstration 

Task 

III 

Program 

Installation 

Task 

IV 

Reports 
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C. Both costs and hours shall he reported for each 
category. Additional reporting categories shall be 
the following: 

Total Contract Costs/Hours 

Fee ■ . , 

Total Contract 

The reporting using the Form 5 33P shall consist of 
columns 7(b) (Planned Value of Work Accomplished) and 
11; (Technical Assessment of Progress). In column 7(b), 
all the above reporting categories shall be included, 
but only dollars are required. The same reporting 
categories shall be used in column 11, except that 
the last two categories above (i.e., Fee and Total 
Contract) shall be omitted. 

D. The monthly reports submission date shall be within 
ten (10) operating days after the closing date of the 
Contractor's accounting month. 

E. Within thirty (30) calendar days after notification of 
approval of the Contractor's proposed final report, 

the Contractor shall distribute a maximum of twenty (20) 
copies of the approved final report as directed in 
writing by the Contracting Officer. In addition, the 
Contractor shall furnish one (1) set of glossy photo- 
graphs of all continuous tone figures included in the 
final report. 

SPECIAL REQUIREMENTS 

Computer Program Minimum Acceptance Standards 

The*, computer program submitted in this contract shall be pre- 
pared for submission to NASA Lewis Research Center in accordance 
with the minimum standards of good practice for computer 
programming as defined below. 

A. The program shall be written entirely in ANSI FORTRAN 
X3.9-1966 and compilable by the UNIVAC FORTRAN V compiler. 
Machine language functions or subroutines shall not be used 
that are not available on the standard FORTRAN V Library. 

B. All subroutine and program options will be fully defined. 

Test input and results will be developed to demonstrate the 
normal program operation. Additional inputs and results 
will be supplied to demonstrate the use of each option. 
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The overall program should be written so as to keep the 
amount of operation card and tape handling to a minimum. 
Operating instructions must be clear. Instructions re- 
quiring halting of the machine, use of external sense 
switches, operator intervention or mounting of tapes after 
program has started execution are not permitted. 

The program description shall include the following: 

1. A complete definition and analysis of the problem 
including equations with engineering units. 

2. A source listing of the program, subroutines, and 
functions, 

3. A complete description of program capabilities, 
logical techniques and options. 

4. A description of all known operational peculiarities, 
limiting cases, and accuracy levels including reasons 
for typical failures. 

5. A list of all input and output parameters with corres 
ponding units and formats. 
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VOLUME I 

COORDINATE FRAME DEFINITIONS 
AND TRANSFORMATIONS 


Consider a cylindrical rolling element bearing with a 

flanged inner ring, Figure Bl. Introduce the coordinate system 

(X,Y,Z) such that the X axis is coincident with the outer 
o o 

ring longitudinal centerline. The Y q axis defines the angle 
^ = 0° and the Z Q axis ^ = 90°. The origin of the (X,Y,Z) 
system is chosen such that the Y - Z Q plane bisects the outer 
ring roller path effective length. The (X,Y,Z) system is fixed 
with respect to the outer ring, and the outer ring center of 
mass is fixed in space. 

Introduce second and third moving coordinate systems, 
(X,Y,Z)j and (x,y,z) c , initially coincident with (X,Y,Z) o but 
attached to the inner ring and cage respectively, and free to 
move through space with these components. 

Introduce a fourth coordinate system (x,y, z)^, at each 
roller location i. This will be referred to as the i-th roller 
inertial coordinate system, (R)^. The origin of the 
system is along the roller axis and in the Y - Z; plane, 

Figure Bl. Note that the x^ axis is parallel with X Q and de- 
fines the roller axis ; y^ is directed radially outward and z^ 
is tangent to the pitch circle. The (x,y,z)^ coordinate frame 
positions are fixed in space, relative to (X,Y, Z) Q . 

Now introduce a fifth coordinate system (5c, y ,.T) - which 
is initially coincident with the (x,y,z)^ but is fixed with 
the roller and free to move through space with it. This frame 
will be referred to as the roller frame (R)^, Figure B2. 
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GENERAL FORM OF THE ORTHAGONAL TRANSFORM 

i It is frequently necessary to express the components of 

a: given vector in at least two coordinate frames. This is par- 

! 

ticularly so in definitions of geometry. Here, initial speci- 
fication is convenient in inertial coordinate frames. However, 
the analysis which follows frequently requires redefinition in 
frames which are in motion. 

The linear orthagonal transformation operator [2], [$TS], 
describes the rotation portion of transformation from the "S" 
coordinate system to the "T" coordinate system. It can be 
written as : 

TS 

[«S] = 0 k/ , 

where k = 1,2,3 and 1 = 1,2,3. 


The most general transformation between Cartesian coordinate 
systems is composed of a rotation and translation. Under trans- 
lation, the coordinates of a point described relative to both 
systems will differ merely by a constant vector, which may be 
defined in either coordinates. The transformation of a vector 
{p}^ whose terminus locates the point P in the "S" coordinate 
frame to a vector {P}^ whose terminus locates the same point 
P in the T coordinate frame can be written as 


{P> T = [ $TS] { {P> s + [$D] {U} E ) S (1) 

where {U}gis the constant vector describing the translation por- 
tion of the transformation, and [$D] is a second rotation trans- 
formation operator. 
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If the vectordJh.in equation (1) is expressed in the final 

iu 

coordinate frame, i.e. E=T, then 

[4>D] = [ $TS ] " 1 (2) 

However, the reciprocal matrix is the transposed matrix 

[$D] = [ 5ts ] (3) 

If translation is specified in the final frame, equations (3) 
and (1) are combined to give: 

{P> T = [ $TS ] {P> s + {U> T (4) 

Conversely, if translation is specified in the initial coordinate 
frame, i.e. E=S, then 

[*D] = [I] (5) 

and 

{p> T = [ns] {p + u> s (6) 

Here, I is the identity matrix. 

Six general transformations need to be related to the inde- 
pendent variables: 

1. Roller inertial (R) coordinate frame to the roller 
(R) coordinate frame, [$RR], 

2. Outer ring (0) coordinate frame to the roller (R) 
coordinate frame, [$R0]. 

3. Inner ring (I) coordinate frame to the outer ring 
(0) coordinate frame, [$0I]. 

4,5,6, Inverse transformations of 1, 2, and 3. 
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j Consider transformation between the roller and roller iner- 
tial frames. Four parameters specify the displacement of the 
roller. They are: two components of translation (u x and Uy-j f 

and j two angles of rotation Yy (roller skew) and y z (roller tilt) 


In terms of the angles of rotation of the i-th roller, y 


yi 


and Y Zi » we obtain: ! 


[ $RR] 


cosY zi sin Yzi 
l-siny . cosy . 

Z ^ U 1 


0 

0 


0 0 1 
which, after expansion, reduces to 


cosYy^ 0 -sinyy^ 

0 1 0 
sinyy^ 0 cosyy^ 


cosY z ^ cosYy^ 
|-sinY z ^ cosYy^ 
sinYy i 


sinY 7 . 

z i 

cos Y 7 . 

z i 

0 


- sinYy^cos Y Z Y 

sinY v . sinY 7 . 
y i z i 


cos Y 


yi 


The inverse transformation is given by 



cosy cos Y v . 

z i >i 

- sin Y Zi 

cos Y v . 
yi 

sinYy i 

= [ $RR] = 

siny_ . 

z i 

cosY, . 

z i 


0 


-slnY y . CO sY z . 

slnY yi 

sinY . 
z i 

cosY yi 


( 8 ) 


Similarly, we obtain these values for the four remaining trans- 
formation operators : j 


[*R 0 ]i - 


[$OR] z = 


1 /| 

0 i 

" j 

°'l 

tm 

1 | 

0 : 


0 

I cos^^ 
j-sin^ 

! 

o 

1 cos^^ 
sin'll 


sin^ 

cos^ 

0 ' 
■sin^^ 
COS^-i 


( 9 ) 


( 10 ) 
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[$0I] - 


[*I0] - 


cosy cosy -siny cosy siny 
z i Yi z i yi 


siny 


cosy„ 

0 

(ii) 


z i 

z i 



"Cosy 
_ z i 

siny 

y i 

siny siny cosy 

z i yi yi j 


c °sy Zi 

cosy yi 

siny . 

z i 

-cosy siny 

z i Ti 


-siny . 
z i 

C0S Yy i 

cosy 

z i 

siny siny 
z i yi 

(12) 

. i 

siny 

yi 

0 

cosy 

yi 





FIGURE Bl: CYLINDRICAL ROLLER BEARING INERTIAL 
COORDINATE FRAMES 
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FIGURE B2 : CYLINDRICAL ROLLER BEARING 
COORDINATE FRAMES. 
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It is essential to the analysis that the geometries of the 
interacting solids are described as accurately as possible. 

This accuracy has to be maintained in parallel \vith description 
generality. 

ROLLER GEOMETRY 

Analysis requires that the roller geometry be represented) 
as a composite of slices which have finite width and radius. 
These radii are computed as functions of roller data. In the 
nomenclature of Figure Cl the k-th roller slice radius is given 
at x = x by 

R k = D eff - C r (1 - coslsin" 1 (£-))) (1) 

in the crowned zone, and by 

R k = D/2 (2) 

in the flat region. Here, 

°EFF = C r + I ' Cc r ' iun) 2 ) h - (3) 

Originally, these radii may be individually specified for a 
completely arbitrary roller geometry. Figure C2 shows eight 
radii which define the geometry of a typical roller. Seven 
"slices" were used in this example. Up to 21 radii and 20 
slices may be specified. 

RING GEOMETRY 

As with the roller geometry, analysis requires that the 
ring geometry be specified in the same "sliced" format. 

In the nomenclature of Figure C3, the k-th ring slice 
radius is given at x = x by - 

R. = i (D' + D + ~) + G m (1 - cos{sin -1 (£ }) (4) 

k,m 2 — — 2 J — m v v G m 
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Here, the upper sign refers to the outer race, the lower to the 
inner race and p^ to the diametral clearance. 

Optionally, these radii may be directly specified for a 
completely arbitrary ring geometry. Figure C4 shows six radii 
which define the geometry of a typical raceway. Five slices 
were used in this example. 

CAGE GEOMETRY 

A minimum of information is required to define the cage 
geometry. The relevant dimensions are shown in Figure C5. 
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FIGURE C3: RING GEOMETRY DETAIL 
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APPENDIX D: LUBRICANT PROPERTY MODELS 


Calculations pertaining to friction on lubricant related 
processes require that the viscosity v and the pressure- viscosity 
coefficient a be known as a function of temperature. 

The kinematic viscosity v (cs) at atmospheric pressure 
is calculated at a given temperature t (°F) from Walther’s 
relation [3] 

log 1(J log 1Q (v + 0.6) = A - B log 1Q (t * 460) (1) 

where A and B are constants determined by evaluating Eq. (1) 
at known values of v = v(t). 


Having calculated v at a specific t, p is computed as 

n ■ vp (2) 


where p the lubricant density at temperature t given by: 

p(t) = p (60 °F) - G { t - 6 0 0 F } (3) 

Here, G is the lubricant coefficient of thermal expansion. 

The pressure- viscosity index, a, is defined implicitly 
by the relation 


n(P»t) = n D e 


ap 


(4) 


However, here a is calculated by a relation developed by 
Fresco [4]: 

(560) 

- 4 ? ( t + 46 0) 2 

a = (2.303) 10 [C+Dlog 10 v+ E(log 10 vrr '(in /lb) 


(5) 


wherein v is evaluated at temperature t (°F) and C, D and E 
are constants tabulated by Fresco. 

Relevant lubricant properties for the oils whose properties 
have been programmed in CYBEAN are listed in Table D1 . 


TABLE D: 
LUBRICANT PROPERTIES i 

Kinematic Walther Equation 

Oil Oil Viscosity (cs) Constants 

No. Type 100°F ^210°F A B 


1 

Mineral Oil 

64 

O 

• 

00 

10.349 

3.673 

2 

MIL-L- 7808G 

12.8 

3.2 

10.215 

3.698 

3 

C- Ether 

25.4 

4.13 

11.452 

4.113 

4 

MIL-L- 2 3699 

28.0 

5.1 

10.207 

3.655 
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VOLUME I 

APPENDIX E: FORCES ACTING ON A ROLLING ELEMENT 

Forces acting on a rolling element can be separated into 
those due to: 

• Elastic Contact - generated at the outer and inner 
ring roller contacts and at the roller end - flange 
contact . 

• Cage Interactions 

• Hydrodynamic Friction - developed in the inlet region 
of the roller/race contact. 

• Elastohydrodynamic - developed within the area of 
rolling-sliding contact of the roller/rings and 
roller/ flanges. 

• Centrifugal Effects 

• Quasidynamic Effects and Element Drag 

The following material develops relationships between forces, 
displacements and speeds for the roller loading described above. 

1. Elastic Contact Forces 

This section deals with the forces that develop as a re- 
sult of the Hertzian elastic interaction between rollers and 
raceways and the roller end and flange. 

la. Roller Raceway Forces 

Roller- raceway elastic forces are a function of the 
relative positions of the rollers and rings. When their 
displacements are specified, interpenetration magnitude 
is computed at various stations along the roller length 
at both the inner and outer rings. These stations are 
actually the disks which would result if the roller- 
raceway effective length were to be sliced by planes 
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parallel to their y-z planes. Slices are assumed to act 
independent of each other. The number of slices and geo- 
metry may either be calculated by the program or input 
as data by the user. (See Volume II and Appendices B 
and C of this volume.) 

The relation between load and deformation for a line 
contact is given by [5], 


10/9 8/9 

Q ■ K* s 


CD 


Here, is a constant, 6 is the elastic deformation and 
is the roller effective length. The k-th slice load 
at the m-th raceway is obtained from 

W. 


«k. 


. r .8/9 10/9 

m m 6k, m 


k 

T~ 

m 


( 2 ) 


The elastic deformation for the k-th slice 6 may 

_ k,nr ' 

be obtained by first computing the vector , which 

connects the k-th roller slice center x-coordinate with 
the k-th ring slice center x-coordinate, Figures E-l and 
E-2. The elastic deformation at the i-th rolling element 
location is then given by 


S kl -l { ®lc ,1 } R 


+ R, - R 


k,l 


+ A 


(3) 


for the outer race contact, and by 


k2 


{B k,2 } R 


‘ R k - 


R 


k , 2 


A 1 

I 


(4) 


for the inner race contact. The terms A 1 and A?: are 

o I 

included in Equations (3) and (4) to account for local 
displacements due to: deviation of the rings from a true 
circular profile and, deflection o.f the outer ring caused 
by interference fit pressure. Particular expressions for 
these effects are presented in Appendix F. 
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Equations (3) and (4) substituted in (2) enable the 
computation of the elastic load magnitude at the k-th slice 
as a function of the independent variables. Direction 
still must be established. 

The unit vectors in the direction of resulting elastic 
force are assumed to act along the lines connecting inter- 
acting slice centers, i.e.: 

jji 

= r x k> m } _ ( 5 -) 

lu kn/R m 1 ,£I , i R J 

l B k,ml 

where C is +1 when m=l and C m = -1 when m=2. 

Equations (2) through (5) are used to give the resul- 
tant force vector acting on each slice of the i-th rolling 


element as 

rryJ D 0 


k,m R 


C { JUS } r 6 ) 

| B k,m| R 


To obtain the vector in the roller frame whose ter- 
minus coordinates will be considered the theoretical point 
of contact, it is assumed that the point lies along the 
line connecting slice centers. The vector 

{P k,m } R = { E k ' C V U k,m } R (?) 

locates the contact point, where g is a vector whose 
terminus coordinates locate the k-th roller slice center. 
Figures El and E2 . 


Using (7) the moment vector generated by elastic 
contact forces is 

{ ^k,m } R = {p k,m } R X { F k,m J R (8;) 
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lb . Roller End- Flange Contact 

Figure E-3 shows a cross section through a cylindrical 
roller bearing having sphere-ended rollers and angled flanges 
In this design, the flange surface is a portion of a cir- 
cular cone whose axis is coincident with the ring x-axis. 
Figure E-4. The apex of the cone is located on the x-axis 
at C. 


The flange surface is described by the equation of 
a right circular cone, 

z = f(x,y) = [Cx - c) 2 cot 2 y f - y 2 ] 1/2 (9) 

where Y_^ is the flange angle, positive for circular cones 
that open in the positive x- direction. 


If we ass tome the elastic contact force generated at 
the roller end flange contact to be normal to the flange 
surface, then its resultant must lie along the line normal 
to the surface z=f(x,y) which passes through the theoretical 

point of contact P(P , P ,P ) . The equation of this line 

X y z 

is given by 


(x-P x ) 

-srr I 


cy-p y ) 

-n~y 


= - (Z 


P ) 


( 10 ) 


9X 'p x ,p y "Sy 'p x ,p y 

Introduce the vector {T)r which locates the roller 
end! sphere origin in the roller frame. We can transform 
the coordinates of this point into the flanged ring frame 
by using the general transform operators listed in Appendix 
B. We shall call the transformed vector {T}j and let 
its terminus point have the coordinates (T x ,t y ,T z ). 


The roller end flange contact elastic force is also 
normal to the roller end surface and its resultant must 
pass through the sphere origin* (T x ,T y ,T z ). Using 
(T x ,T y ,T z ) in equations (9) and (10) we obtain 
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( 11 ) 


(T - P ) = 
v x x J 


' ( T z - p ,) (P x -C) cot- y £ 

[(P x -C) Z cot 2 Y £ -Py] 1/2 


(T - P ) 

y y J 


- P z) P y 

[(P x - C ) 2 COt 2 Y £ -Py] l/Z 


( 12 ) 


P, = [(P v “ C) 2 cot 2 Yf - Pj] 1/2 


x 


(13) 


Equations (11) through (13) contain three unknowns 

(P ,P.P ) and are sufficient to determine the theoretical 
x y z 

point of contact between the roller end and flange. However, 

bjy introducing a fourth equation and unknown, namely, the 

length of the line from points (T ,T ,T ) to (P ,P ,P ) , 

x y z x y z 

we obtain the added benefit of closed form solution. The 

length of a line normal to the undistorted flange surface 

at the point (P ,P ,P ) , which joins this point with the 
x y z 

sphere origin (T x ,T y ,T z ) , is given by 


D - CC T x ' P x> 


(T 


p ) 2 + 

y J 


(T, 


P ) 2 ] 1 / 2 
z J J 


(14) 


After algebraic reduction, the value D is obtained 
from the positive root of the quadratic equation 


n _ -S + (S 2 - 4RT) 1/2 

rjf : : IW :■ -• 


(15) 


where values for S, R and T are 
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R = tan Y £ - 1 


2 sin" Y, 


S = 


[ (T - C) - tanY, (T 2 + T 2 ) 1/2 ] 


cos Y £ lv ‘x £ '-y z 

T = [(T v - C) - tanY, (T v + T 2 ) 1/2 ] 2 


x 


The coordinates P(P ,P ,P ) are given by the following 

a y 4 

closed form functions of D: 


p = 

X 

T y tanY £ 

T z 2 h 

[i + (V> ] t 1 - 
y 

DsinY 

(T y 2 +T 




Dsin y £ 


p y = 

T y 

[1 - 

(T 2 +T z ) 1 / 2 ^ 
v y z 





Dsin y £ 


P z = 

T z 

tl - 

r^TTi/zl 

y z . 



The 

res 

ulting force acts 

along 


7v% 


] + c (16) 


(17) 


(18) 


the flange surface and passes through the sphere center* 
We may use equations (16) through (18) to define the unit 
vector in the direction of this force as 


{ Vl(o) 


where 


cT x 
e 

Id 1 1 

1 e 1 


(19) 


I(o) 


(d M 


I(o) 


!C T x - Vi 
1 V " Vi 
V CT z ' V k ) 


( 20 ) 
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These quantities are used to formulate the expression 
for the forces acting on the e-th flange from its elastic 
contact with the i-th rolling element, 



} 

I 


M’ 



I 


( 21 ) 


where the computation of magnitude, M* , is shown in 
Appendix H. Thtf moment generated by this contact is 
given by 

{ m 1 ) = (PM X { l 1 > (22 

e I ! e i e i ^ 

where 

(P 1 
e 

2 . Cage Analysis 

Several authors [6,7,8] have shown., by experiment and 
analysis, that the cage imposes a substantial effect upon 
rolling element kinematics. Its presence within the bearing 
complement can be especially influential when the bearing is 
operated under light radial load and high speed. Under such 
conditions cylindrical roller bearings can sustain severe skidding 
at the inner raceway contacts. Bearing speeds at any instant in 
time are no longer fully determined by the raceway contact 
forces alone. 

The typical high speed cylindrical roller bearing cage 
is manufactured so that it is supported within the assembly by 
guide rings on its outboard sides. The cage guide rings inter- 
act with shoulders manufactured on either the outer or inner rings. 


} = 
I 


P i 
x ~ 

P j 

y ; 

P k 
z 


(23) 
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The inner ring guide design is prevalent, since the torque induced 
at the guide ring/shoulder interaction tends to be in the direc- 
tion of roller motion and helps in the reduction of skidding. 


The clearance between the guide rails and the ring shoulder 
is sufficient to permit limited extent but full six degree of 
freedom cage motion. However, for current practical considerations 
the cage motion can be considered to reside in a plane normal 
to the bearing axis. This planar motion is encouraged by the 
rotation induced gyroscopic counter moments which tend to keep 
the cage rotation about a single axis. 

... ;• ■ ' . I . 

To analyze the configuration, a coordinate frame (x,y,z) 
is attached to the cage and defined to be coincident with the 
respective ring coordinate frame when the cage is centered on 
the guide shoulders. 


The lubricant forces which develop at the guide ring when 
the cage center is displaced radially to the point y', z ' in 
the y £ rz ( , plane can be conveniently computed from the hydrodynamic 
solution for short, self-acting journal bearings [9]: 


W = + 

y - 


n 0 UL 


Ci-e 2 ) 2 


(24) 



n o lIL 


ire 

4 (1- e 2 ) 3/2 


(25) 


The upper sign is used for inner ring shoulder guided cages 
and the lower sign for the corresponding outer ring condition. 
The value £ is related to the displacement by 


e = (y' 2 + z’ 2 )Vc 


(26) 
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The torque induced on the cage in the presence of lubrication is 
n VR 2 L 

M* = + -°--- „ (2 7) 


Mx* 



— — 

y 


. w z 

d 



1 

0 

0 


M* 

0 

cos^' 

-sin'f'* 

• ■■ ■■ 

W’ y 

0 

sin’l M 

cos^' 

■ 


. W ’z. 


Cl-e 2 ) 1/2 


The forces W'y and W' z act at the point of closest approach 
between the cage guide rail and support shoulder. 

We introduce the azimuthal angle 'l M , which locates the 
point of closest approach and is measured CCW positive from 
the cage ^,-axis. The forces acting on the cage from its inter- 
action with the pilot surface may then be brought into the cage 
reference frame by performing the transformation 


(28) 


The forces which develop within the boundaries of the cage 
pocket are induced by the interaction of the roller with the 
leading and trailing pocket webs and the pocket sides. Magni- 
tudes of these forces depend upon the characteristics of the 
interacting geometries. 

Several pocket designs are currently manufactured for use 
in high speed cylindrical roller bearings. Differences reside 
in the use of materials and the degree of web conformity to 
roller shape. 

Several computation procedures, and thus several opportu- 
nities, exist for the calculation of roller-web interaction 
forces. In the current analysis, the web geometry is simplified 
to that of a rectangular cavity, and the roller to that of a 
cylinder. End effects, such as those which would appear as 
a result of the roller interacting with the cage pocket sides* 
are neglected. 
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j The web forces can be considered to be a function of roller 
rotational speed and the relative position of the roller within 
the cage pocket. Assuming that the rolling element orbital 
velocity remains constant as it traverses the distance 
corresponding to one half the pitch spacing on either side of 
its nominal azimuth position, the value of the roller center 
cage pocket center offset for the i-th rolling element (Z^) 
can be computed from knowledge of the orbital speeds as follows: 

(w ) + (w ) 

j 1 ■ v o . , v o . 

Z. = — ( E [ ! - 1 ]) -z + y' sinijj . - V cosi|k (29) 

l n . 2W i r i c y i v J 

r j = 2 c 

i-2 , 3, . . .n r 

where is the initial circumferential offset and W c is the 
cage rotational speed, assumed to be the average of the roller 
orbital speeds. 

The means for calculating the cage web force evolves from 
the assumption that the pocket load is borne hydro dynamically 
by a lubricant film existing in the gap between the roller 
surface and web. 

The equation relating the film thickness to load magnitude 
is of the form 

F = K'/h o (30) 

where K' is a constant and h Q is the hydrodynamic film thickness. 
In the context of the current numerical analysis, field equations 
containing relationships of this form are extremely difficult 
to solve due to the coupling with many of the independent 
variables. 

A force-displacement relationship of the form 



( 31 ) 
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has currently been chosen. Here, the value of K is selected 
so that equations (31) and (30) are equal at F = 66N (15 lb). 

The latter value is chosen to represent practical upper limits 
for web loading values. Note that the direction of F^ is 
opposite to Z^. 

The radially directed friction force is assumed to be of 
the form 

f x - V F i w 

where is a hydrodynamic friction coefficient and may be com- 
puted from established equations describing the lubrication 
of a rigid cylinder near a plane [10]. The moment about the x-axis 
is given by 

- I \ l F i> < 33 > 

In vector format, the i-th rolling element loading due 
to cage interaction is 


<V >R " 


{He 1 - 
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3. Inlet Region Hydrodynamic Friction Forces 

During high speed lubricated operation, the roller-race 
anjd roller end- flange contacting surfaces are separated by a 
fluid film. The formation of this film depends largely on the 
ability of the contacting bodies to pump lubricant into the 
high pressure region of the contact. 

The fluid immediately preceding the contact, commonly 
referred to as a meniscus, forms a reservoir from which lubri- 
cant is fed. The lubricant is pumped into the high pressure 

i ! i ■■■■■•' 

area by the shearing action created by the rolling/sliding 
motion of the contacting bodies. The non-uniform pressure dis- 
tribution created on the reservoir fluid contributes to rolling 
resistance. 

Expressions for these forces have been developed for point 
(ellipsoidal) contacts [11 , 12 ] , as well as for line contacts [13] , 

The hydrodynamic forces F r , F g and F^ are displayed in 
Figure E5. The special case of contact between two disks is 
illustrated. 

The force F^ acts in the same direction on both contacting 
bodies and opposite to the direction of motion. F acts in oppo- 
site directions on the two bodies in such a way as to tend to 
increase the speed of the slower body and decrease the speed 
of the faster body. 

The forces F 1 and F„ are the integrals of the hydrodynamic 
pressure distribution in the inlet and act through the centers, 
of the two bodies. The component of these forces in the y 
direction represents the (small) portion of the total load sup- 
ported hydro dynamically. The components F nl and F n2 acting 
in the x-direction contribute to the force balance in the 
rolling direction. 
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The magnitudes of F , F s and F n have been obtained in the 
literature [13] . Normal friction forces acting through the cen- 
ters of the two cylinders of radii Ry^ and and of length SL 
are : 


i P x0 C^R*) [h* (R yl *R y2 )]' 1 6 


V 


R y2 

F n2 

■ '■ 


W 


Pumping forces are 


F r = °' 5 < F nl * F n2> 

Sliding friction has been determined to be 

= B 


F s (outer ring) 


F s (inner ring) 


V ( 0uter ) 
'•inner-' 


u ( outer ) 
v inner 7 


- ,outer., 

R r 2nu ( inner )a n 2/3 
2a L K J 


where u = entrainment velocity = (u 2 + u^/2 and V = sliding 
velocity ~ u 2 ~ U 1 

and the dimensionless coefficient P xQ is given in reference [13] 

Retaining vector format, the net force acting on the i-th 
rolling element due to hydrodynamic friction is given by 


(H 1 } 


R 


C p s + F nVl 
( p s + Vy} 
C p s + Vz k 


(36) 


R 


and the net moment acting on the i-th rolling element due to 
hydrodynamic friction is given by 


(F + F )(R ,) i 
{g 1 ) R= |(P r * FslCRylV; 


(37) 


L CF r + F s ) (Ryl ) z kJ _ 
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4. E 1 as tohydro dynamic Forces 

The calculation of forces which develop within the area 
of the concentrated contact follows. 

Pure rolling motion does not exist between the roller and 
raceways. Interacting geometries and elastic contact deformation 
prohibit this. As a result, various points within the contact 
boundaries experience different degrees of sliding. This sliding 
gives rise to fluid shear and asperity interactions, which in 
turn produce forces in the contact area. 

Several authors have addressed this problem during the 
past decade, and produced important empirical and semi-empirical 
models which relate the magnitude of force to measurable variables 

The model employed in the current effort was developed 
by Allen, et al [ 14 ]. 

The traction force is calculated by integrating the fluid 
shear over the entire contact area. Allen assumes the fluid 
shear to behave exponentially with pressure until some critical 
contact pressure is achieved, usually 6.89 N/MM (1,000 psi) , 
above which point the shear can become a constant fraction of the 
normal pressure. This relationship is stated by 


t = n Q e 


ap {?} 


-1 


when t < x 


(38) 


or 


t = Tp 


when n o e ap {s} h -1 > fp (39) 

and >x 


where h is the elastohydrodynamic film thickness, computed 
by the Loewenthal model [15]. 
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Roller raceway contact friction is obtained as the sum 
of the friction forces evaluated at each elemental slice area. 
For a typical slice interaction this is expressed as 


T = 2n 0 h 


e a P H-Cf) 2 ) ■ h’ 1 


where p is the maximum Hertzian contact pressure, b the contact 
half width and z a dummy coordinate measured from the center 
o£ the contact (z = 0) outward. Note that sliding is assumed 
constant across the slice contact area and is computed by simply 
taking the difference between the surface velocities of mating 
slices at the theoretical point of contact (see equation 7) : 


{I* } = {w_} w X { PM 


m R 


km R 


{w o } R X {P km } R 


n X {P* } w 


(41) 


r R 


km R 


Hence, the k-th slice traction force is denoted by 
{P km } R- 

If it is assumed that the resultant slice traction force 
acts through the theoretical contact point, the net moment 
produced is given by 

{U L > " < F km 1 X {T L> R 

The roller end- flange contact is treated in a similar 
fashion, although in this instance the simplifying assumptions 

of uniform sliding/uniform pressure over the contact area is 

» 

invoked, yielding 


{t e } I(o) {s e } I(o) 


t 1 A 1 
e e 


57 


Va***.-**;) 


AL78P022 


Here, is the contact area of the e-th flange contact for 

the i-th rolling element. Sliding velocity is defined as 

the difference between the flange surface velocity and 

roller end surface velocity at the contact point (see equation 23) 




{ 






I + <"o J I 


{Pi>l 


- {w;J; > T x (pM t 
ml el 


(44) 


Finally, the moment produced at the roller end flange 
contact is obtained from equation (23) as 


{u e>I(o) = {F e>I 


{t e } I (o) 


(45) 


5. 


Other Forces 


Rolling Element Drag 

After accounting for all rolling element contact fric- 
tion forces, there is a residual drag acting on the rolling 
element. This force arises as a result of the rolling 
element plowing through the air-oil mixture in the bearing 
cavity. To evaluate this force we employ a model of the 
following form: 

,2 

i -> * 

(46) 


{F DRAG } R " " { t 


PACyd 


8 g 


e i-v 2 

(wj 


] k }1 


where C y is a drag coefficient given in [16). The density 
of the air-oil mixture in the bearing cavity, p, is given 
by 


P = 


p o x 


(47) 


where x is the fractional amount of lubricant in the bearing 
cavity and p Q is the lubricant density. 
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Equivalent inertial forces and moments are introduced 
to apply the equations of static equilibrium when accelera 
tion is present. Using Newton’s second law of motion 
for a single mass, Kellstrom [17] derived expressions 
for these forces and moments. They are: 


i. 



and 




K ) 2 

C"o )2 


dw 1 

o 

dip 



(w 1 ) m D 2 
_o_ 

g 




(48) 


(49) 



t| . 

■ il 
■';$ • 
** -jj 


Note that terms containing mx, my, y s and x are neglected. ** 

Derivatives are approximated from values of speed 
for adjacent elements. 





FIGURE E-l: INNER RING AND ROLLER SLICE 
INTERACTION DETAIL. 
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FIGURE E-4: ROLLER END-FLAI 
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NOMENCLATURE 
Appendix E 

? 

A frontal area of a rolling element {in ) 

C cage land radial clearance (in) 

cage pocket radial clearance (in) 

D roller diameter (in) 

E modulus of elasticity (PSI) 

I effective cross section moment of inertia (in^) 

L cage land width (in) 

load deformation constant for a line contact 

K load deformation constant for a point contact 

p 

R a/R > . t + 1/Ry 2 ) (in) 

R c cage land radius (in) 

R^ radius of the k-th roller slice (in) 

R, radius of the k-th raceway slice of the 

,m m-th raceway (in) 

U entrainment velocity (in/sec) 

V net sliding (in/sec) 

W slice width (in) 

b contact half width (in) 

d pitch diameter (in) 

g gravity constant 

h film thickness (in) 

l roller raceway effective length at the m-th 

m raceway (in) 
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m 

roller mass (lb-sec /in) 


n i 

number of roller-raceway slices 


n r 

number of rollers 


r s 

roller end sphere radius (in) 


u 

surface velocity (in/sec) 


w 
i m 

rotational speed of the m-th raceway 

(rad/sec) 

w 1 

0 

orbital speed of the i-th roller (rad/sec) 

< 

rotational speed of the i-th roller 

(rad/sec) 

H 

azimuth angle locating the i-th rolling element 
measured CCW positive from the y-axis in the y- 
plane of the outer ring frame 

Y f 

flange angle (rad) 


Y y 

roller skew (rad) 


Y z 

roller tilt (rad) 


% 

absolute viscosity (lb-sec/in ) 


{V} 

a 

a vector in the a coordinate system 


A 

iV >a 

a unit vector in the a coordinate system 

imj 

magnitude of the vector V 



i,j,k unit vectors in the x, y and z directions, respectively 

Program input and output are in SI units although calculations 
are performed in English units. 

Subscripts 


e 

refers to the e-th flange 


i 

refers to t.he i-th roller in complement 


k 

refers to the k-th slice 


Z 

refers to the il-th rolling element outer 
contact load 

ring 

m 

refers to the m-th raceway, i.e., m=l for 
outer race, m=2 for the inner 

the 

Notation 



X 

vector cross product 
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VOLUME I 

.APPENDIX F; FITTING OF NON-CIRCULAR RINGS — 

INTO NON-CIRCULAR HOUSINGS 

| Several authors have shown that the flexural stiffness 
of the outer and inner rings imposes a first order effect on 
the rolling element load distribution, and hence on the bearing 
performance. However, these investigations have been limited to 
the analysis of bearings for which the reactive housing/ring 
load distribution has been assumed. In the current analysis the 
reactive load distribution and hence the deformed shapes of 
both the housing and outer ring are left as variables to be 
defined by solution. 

The rolling element load at the outer ring is substantially 
higher than that present at the inner ring for high speed bearing 
applications. This is due to roller forces generated by centri- 
fugal effects. High speed cylindrical roller bearings in parti- 
cular, are typically operated under light imposed radial loads. 
Thus, for the present analysis, the deformation of the composite 
outer ring/housing is detailed and the inner ring is assumed to 
be a flexurally rigid elastic solid. 

To induce roller preload, many cylindrical roller bearings 
are manufactured with a non-circular outer ring surface.. 

The radius of the outer profile of such a ring is a 
function of azimuthal position, i.e.: 

rj (0) = Rj (1 + e j cos n j (0 - $j) ) (1) 

Additionally, the housing itself may also be manufactured in 
a similar manner so that: 

r (0) = R (1 + e cos n_ (9 " $ )) (2) 

o ^ o v o o v o'' - 

Here, R Q and Rj are the average radii of the housing and ring, 
e Q and ej the non-circular semi -amplitude (expressed as a frac- 
tion of the average radius) , and are phase angles and n Q 
and n^ the number of lobes (Figure F-l). 
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The deformed shape that the outer ring assumes is governed 
by ; the free shapes of the ring and housing, the rolling element 
loads and the reaction pressure distribution at the housing/ring 
interface . 

Currently, no closed form solution exists for the defor- 
mation of thin circular rings loaded by an arbitrary pressure 
distribution. However, expressions have been derived for the 
deformation of a ring under the influence of discrete point 
radial loads. These solutions may be used to construct an equi- 
valent load distribution on a given ring. The deflection of 
a particular point subject to an arbitrary load distribution 
is then defined by the superposition of the deflections obtained 
from each point load in the distribution. 

Additionally, from equilibrium considerations, we require 
that the sum of the reaction loads has to equal the sum of the 
roller applied loads: 

m n 

Z F.cos0. = Z G-cosy. (3) 

i=l 1 1 j=l 3 3 

m n r 

Z F.sin0. = Z G.siny. (4) 

i=l 1 1 j-1 3 _ 3 

Where F^ , i=l, m are the roller loads and Gj , j«l, n are the 
point loads which approximate the reaction pressure distribution. 

The deflection of the ring from both the roller and reac- 
tion loads can be computed from: 

m n 

D t (0) = Z F.C t (B-0.) - Z G.C (e-y,) (5) 

1 i=l 1 1 1 j=l 3 1 3 

where C (0) is defined as the outward deformation of the ring 
at 0=3 due to an outward unit load at 0=0. 
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Here, the terms C-j-(g) are referred to as influence coeffi- 
cients of the ring and their explicit form may be found in the 
published literature [18]. 


Similarly, the effect of reaction loading upon housing 
deformation is computed as 


D 0 (3) = 


n 

£ 


j-1 


G jC 0 (e-Yj) 


C6) 


where C Q C 3 D carries an analogous definition to Cj(g). Note 
that the influence coefficients for a housing of completely 
arbitrary geometry can be obtained by the use of a finite ele- 
ment analysis. 

Additional constraints are needed now to confine the de- 
formed ring shape within the physical boundaries of the deformed 
housing. Having defined the set of n angles of interest, 0j(j = l,n), 
the problem is reduced to one in which we must find a subset 
of indices, (k=l,n') such that 

at 0 k = 0 £ : D q - Dj + A = 0; k = l,n’ (7) 

k 

and at all other angles, 0^=0^: 

D 0 - Dj + A > 0; j = l,n (8) 

- 1 * \ 

■ ■ i 

where A is the clearance between the undeformed ring outer sur- 
face, whose center is displaced to coordinates (Xj, y^.) , and 
the housing: 

A = R q (1 + e Q cosn 0 (0 - ; $ Q )) - Rj (1 + ej-cosnj (0 - $j)) 

x (9) 

+ Xj cos0 + yj sin© 


70 


AL78P022 


Equations (5) through (9) are combined and evaluated at 
anj angle ^ to yield 

n m n 

D -D t 46= 2 G.C (0-0. ) - 2 F. C T ($-0 - ) + 2 G,C T (8-9 ) 

0 1 j = l J ° \ i=l 1 1 1 j=l J 1 

■ ■ I ' ' ■ . - . . . . ; ' ... • ■ - 

+ R (1 + e cos (n 8) ) - R T (1 + e T cos (n T 6)) 


( 10 ) 


+ XjCosB + yjSin8 


Whbn a set of (k = l,n') is assumed, equation (10) evaluated 
atj 3 = ®£. (i=l,n’) and equated to zero and equations (3) and 
(4) witli form a set of n'+2 linear equations in the 

n'+2 unknowns G,(k=l,n'), Xj and y^. After solution of that 
set of equations, we must check that all G^ (k=l,n') are positive 
then equation (8) must be checked at all 0^ with j=l,n and with 
(k=l,n') and must be positive. 


i - ■ 


1 

( 

i 
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VOLUME I 

APPENDIX G: GOVERNING EQUATION SET AND SOLUTION ALGORITHM 


EQUATIONS OF EQUILIBRIUM 

The determination of roller, ring and cage loading detail 
during operation, requires that the representation of these com- 
ponents satisfy equations of force and moment equilibrium. A 
total of 6n„ + 5 equations may be formulated. 


Jjjj. 

I 

•*** ft . 

’ ’i '■■it 
■ s*4 K . 

ii-'i. I 


The formulation of the equations of equilibrium in turn, 
frequently requires the specification of a force or moment in 
a coordinate frame different from the one in which it was ori- 
ginally determined. This is accomplished by vector manipulation 
and the use of the transformation operators described in 
Appendix B. 

INNER RING EQUILIBRIUM EQUATIONS 

The computation of roller contact elastic and friction 
forces enables the writing of inner ring equilibrium .equations 
by the summation of these loads over all rolling element 
locations i : 




•js 


n 


r 

I 


n 

[ 2 


i=l k=l 


c {F i,2 +T i,2 v 


+ ♦ < p Ah = ° 

e = l 


(1) 


t 


j, I.V'4,2 * U k,2>d * + ;M A>I 

1=1 lc-1 ’ ’ e=l 


= 0 


( 2 ) 


Here, {F A > and {M^} represent the specified applied loading. 
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; Hydrodynamic loading is excluded from this equation set 
because of its extremely small magnitude when compared with con 
tact and applied loads. The six vector component equations (.1) 
and (2) with the appropriate imposed boundary values define the 
inner ring equilibrium field equation set. However, in bearing 
design analysis it is very convenient to specify the radial 

components of {F.} and a misalignment angle. In this case the 

■**' I 

moments decouple from Equation (2) and consequently may 

be evaluated directly. The remaining equations of equilibrium 
are reduced to the Yj and Zj components of (1). 

ROLLING ELEMENT EQUILIBRIUM EQUATIONS 

The (6n r ) rolling element equations of equilibrium are 
obtained by summing all forces and moments applied to the i-th 
roller as follows : 


2 n 


V ■ E , {F k,m + T lU¥ * 

m=l k=l ’ e-l 


{F 1 } + {H 1 } + {Fj.} + { p n rag }_ = 0 

r q R Urag R 


C R 



m ; 


2 n . 2 . 

£ ( E m + U V m } ) + Z ) + 

= 1 k=l K ’ m k ’ m R e-l e R 


(4) 


{M 1 } + {G‘ l } + {MM = 0 


r 


ui 


C R 


R q R 
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CAGE EQUILIBRIUM EQUATIONS 

The cage maintains equilibrium by balancing its pocket and 
land forces : 

n r 

Z [ (-f - 3 sin#. + (- F. ) cos#.] - W = 0 (5) 

i=l 1 . y 

n r 

Z [(-£•) cos#. + (F. ) sin#. ] - W =0 (6) 

i=l 1 11 i z 


X X 

M + i Z C-f.) - Z M* =0 (7) 

x z i=l 1 i=l 

SOLUTION ALGORITHM 

Several computation procedures, and hence several oppor- 
tunities, exist for the economic calculation of values for the 
set of independent variables. In contrast to the generality 
maintained in the bearing analysis formulation, specific solu- 
tion algorithms are required to address the specific problems 
arising during the computation of bearing performance detail. 

The following' physical conditions make the field equation set 
particularly difficult to solve: 

1) Con tact/No- Contact - The rollers transmit force 

through physical contact with the outer and inner 
rings and one or more flanges. It is known that the 
rollers must at high speeds, contact at least the 
outer ring at all times (centrifugal force) «, However, 
similar knowledge of flange contact kinematics is 
unavailable, thus giving rise to very poor initial 
guesses of skew and axial position. 
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2) Piece-wise Continuity of Function Derivations - As a 
result of clearances in the bearing, such as roller 
end play and radial looseness, the governing equation 
set becomes piece-wise continuous in its derivatives. 
This is especially true of the inner ring, roller 
radial, skew moment and axial equilibrium equations. 

3) Accuracy - Loss of accuracy, through round-off errors, 
is especially influential when computing approximations 
to the derivatives of the field equation set. 

A solution algorithm was constructed which minimizes these 
problems. Equilibrium is sought through automated Newton- Raphson 
iteration. Automatic variable damping is included to eliminate 
oscillations in localized minima of the function hyperspace. 

An automated field equation set partitioning scheme was 
developed so that many different subsets of equations and inde- 
pendent variables could be experimented with to determine an 
optimum solution algorithm. Experimentation with this parti- 
tioning scheme led to the isolation of axial displacement and 
roller skew as the most difficult variables to control, along 
with the axial and skew moment equilibrium equations, the most 
difficult to equilibrate. Field equations written in the stan- 
dard "force-displacement" format could not be brought to numeri- 
cal equilibrium. 

This was caused by the high sensitivity of the flange in- 
duced force to variation in axial and skew positions. Large 
roller end radii need only very small penetration into the 
flange to yield high reactive forces. 
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To avoid this numerical problem, an alternate solution 
algorithm was devised. Raceway induced roller loading is com- 
puted in the usual manner using standard "force- displacement" 
relationships. The net axial raceway induced load, F^, required 
to produce equilibrium is obtained and specified to equal the 
axial component of the flange induced load. The contact side, right 
or left flange, is obtained by noting the sign of F NX » Using 
the unit vector in the direction of the resulting flange force, 
which was given in Appendix E, along with the EHD traction model, 
a flange induced load vector and contact location is obtained. 
Moments are computed using the vector cross product. 

The above calculated flange force vectors, together with 
the raceway induced forces, provide roller equilibrium in the 
axial direction. With this equilibrium obtained, the axial 
equilibrium equation is removed from the field equation set and 
replaced by the equivalent kinematic relationship which states 
that the flange penetration is zero (see Equation (15) , Appen- 
dix E) i . e . : 

< b 1 = D ± - r s = 0 (8) 

Here, r is the roller end sphere radius and the numeric sub- 
script refers to the specific side on which contact occurs. 

Using the above algorithm, one establishes contact at one 
flange and satisfies axial equilibrium. 

This algorithm was subsequently expanded to include two 
contacts. By comparison of the flange induced roller moment 
with raceway induced moments, one may distribute the flange 
induced loading between two flange contacts (if necessary) such 
that both axial and skew moment equilibrium is achieved. 
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Consequently, the skew moment equilibrium equation is dropped 
from the equation set and replaced by the equivalent kinematic 
equation : 



Full algorithm detail and logical options are shown in Figure G-l 
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NOMENCLATURE 
APPE NDIX G; 

elastic roller/raceway contact load (lbs -see equation 
6, Appendix E) 

friction force at roller/raceway contact (Ibs-see 
equation 40, Appendix E) 

roller end/flange contact load (lbs) 

elastic roller/raceway contact moment (in-lbs 
see equation 8, Appendix E) 

friction moment at roller/ raceway contact (in-lbs 
see equation 42, Appendix E) 

m* friction moment at roller end/flange contact (in -lbs) 

F 1 roller/cage web contact load (lbs - see equation 34, 

c Appendix E) 

H 1 hydrodynamic roller/raceway contact force 

(lbs - see equation 36, Appendix E) 


*k,m 


k,m 

i 

e. | 

i 

k,m 


U 


k,m 


quasidynamic roller force (lbs - see equation 48, 
Appendix E) 


F 


l 

Drag 


roller drag force (lbs - see equation 46, Appendix E) 


M 1 moment generated by roller/cage web contact 

c (in-lbs - see equation 35, Appendix E) 


G 1 moment generated on roller by hydrodynamic forces 

(in-lbs - see equation 37, Appendix E) 





M. 


x * 


quasidynamic roller moment (in-lbs - see equation 49, 
Appendix E) 

radially directed cage pocket/roller contact friction, 
force (lbs - see equation 32, Appendix E) 

cage pocket/roller contact load (lbs - see equation 
31 , Appendix E) ■/ \ 

cage/pilot surface forces (Wy, W z -lbs) and moment 
(M x -in-lbs) about cage axis of rotation (see equation 
28, Appendix E) 
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{V) a a vector in the a coordinate system (coordinate frames 

u are defined in Appendix B) 

Subscripts /Super scripts 

e refers to e-th flange, e = 1,2 

i refers to i-th roller, i = 1 , 2 , 3, . . . , 50 

k refers to k-th slice (per Appendix E - Roller Raceway 

Forces) k = 1,2,3,. ..,20 

m refers to m-th raceway, m=l for outer ring, m=2 for 

inner ring 
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VOLUME I 

APPENDIX H: BEARING FATIGUE LIFE CALCULATIONS 


INTRODUCTION 

Within CYBEAN, roller bearing fatigue life is calculated 
using Lundberg-Palmgren [ 13, 20] methods. The value computed is 
then modified by multiplicative factors which account for 
material and lubrication effects. 


ROLLER BEARING RACEWAY LIFE 

The Load distribution across a line contact is represented 
by a number of slices. The L^q fatigue life of a given slice 
is : 


J 10mk 


r - ( ^cmk -' 4 
Qemk 


CD 


Here, Q cm j c is defined as the load for which a slice will have 
90 percent assurance of surviving 1 million revolutions. Letting 
m refer to a raceway, k to a slice, where n is the index of the 
last slice, the explicit form of Q cm ] c is given by: 


49500 A { D ( 1 ± y)} 1 * 074 £ 0,778 ( Y j0.222 

m m 


•cmk 


( N(1 + y n )} 


0.25 


( 2 ) 


Here, N is the number of rollers and, the upper sign is used for 
the outer race, the lower sign for the inner race. 

Q em k is the equivalent load for the slice: 


n 


Qemk n ^ ^mkj ^ 
r j =1 J 


1/e 


CD 
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Q m kj is the individual roller contact load or the k-th slice and 


e = 4.0 or e = 4 . 5 , depending respectively upon whether the 
applied load rotates or is stationary with respect to the race- 
way in question. 


L 1q life of a raceway is giv- n by 


n 


L 10m " a 2 a 3 a 3 a 4 ! L;LOmk ^ ^ 


(4) 


where e is the Weibull slope exponent, here taken to be 9/8 


is a life improvement factor to account for improved 
materials 


a^ is a life improvement factor to account for full film 
lubrication 

* 

a^ is less than 1 when full film lubrication is not 
obtained [21] . 

a^ is a factor which accounts for materials having a 

modulus of elasticity other than that of basic steel. 


BEARING LIFE 

L, n life of the bearing considering both raceways is: 
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NOMENCLATURE 
APPENDIX H , 

' »/ d n, 

D roller diameter (in) 

d pitch diameter fin) 

m 

l slice width (in) 

X capacity reduction factor 

X^ = •& when K=2 , 3, . . . , (N-l) 

X^ = -53 when K=1 or N, N=Number of slices 

number of rollers 
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VOLUME I 

APPENDIX I: THERMAL MODEL 


TEMPERATURE CALCULATIONS 

After each calculation of bearing generated heat rates , 
either steady state or time transient temperature analysis may 
be performed. The computations are terminated in the following 
manne r : 

1. The steady state case terminates when each system 
node temperature is within A “Centigrade of its 
previously predicted value,, The value for A is 
specified by the user. 

2. The transient calculation terminates when the 
user specified computation time is reached or 
when one of the system temperature nodes exceeds 
600 °C (1112 °F), 

STEADY STATE TEMPERATURE MAP 

The physical structure is considered to be divided into a 
number of elements represented by nodes. Heat flow to node i 
from surrounding nodes j, plus the heat generated at node i, 
must equal zero to satisfy the definition of steady state 
conditions . 

After each calculation of bearing generated heat, which 
results from a solution of the bearing portion of the program, 
a set of system temperatures is determined which satisfy: 

q. = q . + q . = 0 for all nodes i (1) 

i oi n gi 

where q . is the heat flow from all neighboring nodes to node i 

q . is the heat generated at node i. Values ere 

calculated to represent heat Created by bearing 
friction. 
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The resulting set of field equations is solved with a modi 
fied Newton-Raphson method which successfully terminates when 
either : 


At, 


— <. 6 for all nodes i 


( 2 ) 


here, At represents the Newton-Raphson correction to the i-th 
tempe rature . 


or; 


i n C EQ ) . \/2 

1 y r L < a 

100 . \ 1 n J “ 

i = l 


(3) 


where, n _ num tj er 0 £ field equations 

EQj = residue of the i-th field equation 

6 = .001 

TRANSIENT TEMPERATURES 


The net heat q^ transferred to the i-th node heats the 


element , i . e . : 


dt . 

pC p, v i ht7 


w 


where p = density 

Cp = specific heat 
V = volume of the element 
t = temperature 
T = time 


The temperatures, t 0 i> at the time of initiation T=T g are 
assumed to be known, that is 


t oi = t. (TJ 
1 v s J 


i 1, 2 , . . . , 
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C l 



fc £ 


The problem of calculating the transient temperature dis- 
tribution in a bearing configuration thus becomes a problem of 
solving a system of nonlinear differential equations of the first 
order with prescribed initial conditions. The equations are 
nonlinear since they contain radiation terms and free convection, 
which are nonlinear with temperature as will be shown later. 

The simplest and most economical way to arrive at a solution is 
to calculate the rate of temperature increase at the time T = T^ 
from equation (4) and then compute the temperatures at time 


T^ + AT from 


dt 

t k+ 1 t k + W 


AT = t 


k PC V 
P 


AT 


( 6 ) 


r 

« * 


CALCULATION OF TRANSIENT TIME STEP 


too large, the tempera- 


If the time step AT used is chosen 
tures will oscillate; if it is chosen too small the calculation 
will be costly. It is therefore desirable to choose the largest 
possible time step that does not give an oscillating solution [22,23] 


dt i,k+l 

~ dt i,k 


1 2 

x > u > • • 


n 


(7) 


If this derivative were negative, the implication would 
be that the local temperature at node i has a negative effect 
on its future value. This would imply that the hotter a region 
is now, the colder it will be after an equal time interval. 
Oscillating solution would result. 

Differentiating (6) for node i and combining with (7), 
the time step size condition is 


dt i,k + l _ , 


AT i 


P,C V. 

1 l 


dq, 


clt7 - 
i 


> 0 


(8) 
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The derivative dq^/dt^ is approximated using the forward 
difference operator 



g i C t i + At i )-q i (t i ) 

t i 


( 9 ) 


The values A T^ which satisfy the equality in Equation (8) 

are calculated. The array is searched, and a value of AT 

rounded off to one significant digit smaller than the smallest 

of the At. obtained is used, 
x 


CALCULATION OF HEAT TRANSFER RATE 


Heat transfer mechanisms which occur in a bearing 
application are: 

• Conduction between inner ring and shaft 
and between outer ring and housing. 

• Convection between the surface of the 
housing and the surrounding air. 

• Radiation between the shaft and the housing. 

• Forced convection between the bearing and 
circulating oil. 

All the above mentioned modes of heat transfer are con- 
sidered in calculations of the heat balance at each given node. 

GENERATED HEAT 

A heat source may exist at node i. The quantity represen- 
ting the source magnitude must be added to the net heat flowing 
from neighboring nodes. 
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When the heat source is other than a cylindrical roller 
bearing, it may be considered to produce known amounts of power, 
in which case constant numbers are entered as input to the 
program (see example in Volume II) . 


CONDUCTION 


The heat flow q • • which is transferred by conduction 

^ 1 9 J 

from node i to node j, is: 


l cx,3 


TT: C'i'V 


( 10 ) 


where X = the thermal conductivity of the medium. 
H = length between i and j 

FREE CONVECTION 


Free convection between a solid medium and air, the heat 

q.vi,j trans£f 

from the equation 


flow q - • transferred between nodes i and j can be calculated 

VI , J 


Hvi.j - v I l r t j ; • SI ™ “i-'j' 


(ii) 


where 


a y = the film coefficient of heat transfer by free 
convection 

A = the surface area of thermal contact between the media 

d = is an exponent, usually = 1.25, but any value can 
be specified as input to the program 


SIGN = X 


1 ti > tj 


;1 if < tj 


The value of can be calculated for various cases [22] 
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FORCED CONVECTION 


Heat flow q . - transferred by forced convection can be 
l wi,j 

obtained from the following equation. 


q . . = a A(t. - tO 
M wi,j w 1 2 


( 12 ) 


where a is the film coefficient of heat transfer during forced 
w 

convection. This value is dependent on the actual 
shape, the surface condition of the body, the difference 
in speed, as well as the properties of the liquid 
or gas. 


In most cases, it is possible to calculate the coefficient 
of forced convection from a general relationship of the form, 


N u ■ aR e P T 


(13) 


where a, b, and c are constants obtained from handbooks, [24], 

R and P are dimensionless numbers defined by 
e r ' 


N u = Nusselt number = 

L = characteristic length 
1 = conductivity of the fluid 

R e = Reynold's number = ULP/n 
U = characteristic speed 
P = density of the fluid 
n = dynamic viscosity of the fluid 
P r = Prandtl's number = PCL /A 
C„ = specific heat 


CYBEAN can accept a specified constant value for the coefficient 
of convection, or, let it vary with actual temperatures. The 
variation is determined by the viscosity variation. Input can 
be given in one of three ways, for each coefficient. 
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Constant Viscosity 

1. Values of the parameters in Equation (13) are given 
as input and a constant value of is calculated 
by the program. 

Temperature Dependent Viscosity 

2. The coefficient a for turbulent flow and heating 
of petroleum oils is given by 

a = k ;, tn (t) ) k 10 (14) 

j w g ^ j J 

where and are given as input together with 
viscosity at two different temperatures. 

3. Values of the parameters in Equation (13) are given 
as input. Viscosity is given at two different 
temperatures . 

RADIATION 

If two flat, parallel surfaces of same surface area A, are 
placed close together, the heat transferred by radiation between 
nodes i and j representing those bodies, will be, 

q Ri , j = e °A [(t i +' 273) 4 - (t. ;+ 2 73) 4 ] (15) 

here, £ is the surface emissivity, and c is the Stefan- 
Boltzmann radiation constant. 

Heat transfer by radiation under other conditions cari also 
be calculated ['22]. The following equation, for instance, applies 
between two concentric cylindrical surfaces: 

eoA. [(t. + 273)* - (t. + 273) 4 ] 
qRi »j ’ 1 + (1-e) (A i /A e ) C16) 


AL78P022 

where is the area of the inner cylindrical surface 

A g is the area of th^ outer cylindrical surface 

FLUID FLOW 

Between nodes established in fluids, heat is transferred 
by transport of the fluid itself and the heat it contains. 





r' / / 

FIGURE 1-1 



/ / / 


FLUID HEAT NODES 





Figure 1-1 shows nodes i and j at the midpoints of conse- 
cutive segments established in a stream of flowing fluid. 


The heat flow q . . through the boundary between nodes 

111 j J 

i and j can be calculated as the sum of the heat flow through 

the middle of the element i, and half the heat flow q . trans- 

l oi 

ferred to node i by other means, eg., convection. 

The heat carried by mass flow is. 


q _ . = p. C V • t • = K • t 
l fi 1 p. 1 1 


l i 


(17) 


where = the volume flow rate through node i. 


The heat input to node i is the sum of the heat generated 
at node i (if any) and the sum over all other nodes of the 
heat transferred to node i by conduction, radiation, free and 
forced convection. 
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q oi q G,i + Z ^ q ci,j + q vi,j + q wi,j + q Ri,j^ 


The heat flow between the nodes of Figure I -1 is then. 


q ui,j q fi + q oi^ 2 


(19) 


If the flow is dividing between node i and j, as in Fig. 1-2, 
then the heat flow is calculated from 


q . . = K. . (q £ . + q . /2) 
n ux,j lj ^ n fi n oi 


( 20 ) 


where = the proportion of the flow at i going to node j , 

0 < K. . <1. 
ij “ 



TOTAL HEAT TRANSFERRED 

The net heat flow rate to node i can be expressed as, 


q i = 


q G,i + 


z (q 


ci,j 


l Ul,J 


+ q . . + q . 

H vi,j h wi,3 


l Ri , j 


i) ( 21 ) 


The summation should include all nodes j , which interact 
with i. Unknown temperatures as well as those specified as 
known should be included. 
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CONDUCTION THROUGH A BEARING 

The conduction between two nodes is governed by the thermal 
conductivity parameter X. The value of X is specified at input. 

An exception occurs when one of the nodes represents a 
bearing ring and the other a set of rolling elements. Here, the 
conduction is calculated separately using the principles described 
below. 

THERMAL RESISTANCE 

It is assumed that the rolling speeds of the rolling elements 
are so high that the bulk temperature of the rolling elements is 
the same at both the inner and outer races, except in a volume 
close to the surface. The resistance to heat flow can then 
be calculated as the sum of the resistance across the surface 
and the resistance of the material close to the surface. 

The resistance ft is defined implicitly by 

At = ft • q (22) 

where At is the temperature difference 

and q is the heat flow 

The resistance due to conduction through the EHD film is 
calculated as 

ft x = (h/X) • A (2 3) 

where h is taken to be the calculated plateau film 
thickness 

A is the Hertzian contact area at the specific 

rolling element-ring contact under consideration. 

X is the conductivity of the oil. 
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The geometry is shown in Figure I- 3(a). Asperity conduction 
is not considered. 

So far, a constant temperature difference between the sur- 
faces has been assumed. But during the time period of contact, 
the difference will decrease because of the finite thermal 
diffusivity of the material near the surface, Figure I- 3(b). 

To points at a distance from the surface, this phenomenon 
will have the same effect as an additional resistance acting 
in series with ft^. 

This resistance was estimated [25] as, 
n 2 

where £ „ 
re 


= 1 f-li-i 1 / 2 C241 

-re.i :l i V 

= contact length, or in the case of an elliptical 
contact area, 0.8 times the major axis 


A = heat conductivity 

\p = thermal diffusivity = A/(p» C ) 

p = density 



specific heat 


b = half the contact width 


V = rolling speed 


The resultant resistance is 


ft = ft, + ft~ 

res 1 2 


(25) 


There is one such resistance at each rolling element. 
They all act in parallel. The resultant resistance, ^ res » is 
thus obtained from 

: ■ i n i 

**• _ r> J- r 


1 g i 

^res i=l ^res,i 
98 


(26) 
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direction of rolling 


(b) Temperature Distribution at Rolling, 
Concentrated Contact Surfaces 

FIGURE 1-3: CONTACT GEOMETRY AND TEMPERATURES 
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ATTACHMENT 1 


Roller Bearing Specifications 


Inner Race 


Bore Dia. 

mm 

(in) 


118 

(4.6457) 

Raceway Dia. 

mm 

(in) 


131.66 

(5.1834) 

Flange Dia. 

mm 

(in) 


137.47 

(5.4122) 

Width 

mm 

(in) 


26 .92 

(1.060) 

Groove Width 

mm 

(in) 


14.59 

(. .5746) 

Flange Angle 




. 6 

deg. 

■ Race 






Outer Dia. 

mm 

(in) 


164.49 

(6.4760) 

Raceway Dia. 

mm 

(in) 


157.08 

(6.1842) 

W i dth 

mm 

(in) 


23.9 

( .942) 

:rs 






Diameter 


mm (in) 


12 .65 

( .4979) 

Length-overall 

mm (in) 

- 

14.56 

( .5733) 

-effective 

mm (in) 


13.04 

( .5133) 

- flat 


mm (in) 


8.40 

( .3307) 

Crown Radius 


mm (in) 


622.3 

(24.5 ) 

End Radius 


mm (in) 


381.0 

(15. ) 

Number 




28 


Land Diameter 



mm (in) 

137.95 

(5.4312) 

Axial Pocket 

Clearance 

mm (in) 

. 02 0 ( .0008) 


Cage 


Tangential Pocket Clearance mm (in) 
Single Ball Width mm (in) 

Operating Conditions 

Shaft Speed 

Bearing Radial Load 4450 N 

Oil Inlet Temperature 366. 5I( 

Misalignment of Races 


. 2 2 1 ( .0087) 
4.6 ( .18 ) 


20,000 rpm 
(1000 lb) 
(200 °F) 


0 deg § .5 deg. 
.75 deg. 
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